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We study the coupling between backward- and forward-propagating wave modes, with the same group ve- 
locity, in a composite right/left-handed nonlinear transmission line. Using an asymptotic multiscale expansion 
technique, we derive a system of two coupled nonlinear Schrodinger equations governing the evolution of the 
envelopes of these modes. We show that this system supports a variety of backward- and forward propagating 
vector solitons, of the bright-bright, bright-dark and dark-bright type. Performing systematic numerical sim- 
ulations in the framework of the original lattice that models the transmission line, we study the propagation 
properties of the derived vector soliton solutions. We show that all types of the predicted solitons exist, but 
differ on their robustness: only bright-bright solitons propagate undistorted for long times, while the other types 
are less robust, featuring shorter lifetimes. In all cases, our analytical predictions are in a very good agreement 
with the results of the simulations, at least up to times of the order of the solitons' lifetimes. 



PACS numbers: 41.20.Jb, 42.65.Tg, 78.20.Ci 

I. INTRODUCTION 

Left-handed (LH) metamaterials are artificial, effectively 
homogeneous structures, featuring negative refractive index 
at specific frequency bands where the effective permittivity e 
and permeability ji are simultaneously negative (THSj. In fact, 
all known realizations of LH metamaterials rely on the use of 
common right-handed (RH) elements and, thus, in a realistic 
situation such a composite material features both a LH and a 
RH behavior, in certain frequency bands. Physically speaking, 
the difference between the two is that in the LH (RH) regime, 
the energy and the wave fronts of the electromagnetic (EM) 
waves propagate in the opposite (same) directions, giving rise 
to backward- (forward-) propagating waves. 

Transmission line (TL) theory constitutes a convenient 
framework for the analysis of LH metamaterials. Such an 
analysis relies on the connection of the EM properties of the 
medium (e and fi) with the electric elements of the TL's unit 
cell, namely the serial and shunt impedance. As mentioned 
above, in practice composite right/left-handed ( CRLH) struc- 
tures are quite relevant, giving rise to pertinent CRLH-TL 
models. These models are, in fact, dynamical lattices which 
can be used for the description of a variety of metamaterials- 
based devices and systems, such as resonators, directional 
couplers, antennas, etc fD^]. 

Nonlinear CRLH-TLs, with a serial or/and shunt 
impedance depending on voltages or currents, have also 
attracted attention. Such structures may be realized by 
inserting diodes - which mimic voltage-controlled nonlinear 
capacitors - into resonant conductive elements (such as 
split-ring resonators) fSllTll. Such nonlinear CRLH-TL 
models have been used in various works dealing, e.g., with 
the parametric shielding of EM fields | 8 |, the long-short wave 
interaction |9|, or soliton formation |10-12|. Experiments 
in nonlinear CRLH-TLs have also been performed (see 
the review LI 3 J ), and formation of bright |il4| il5J or dark 



(TsKISl envelope solitons, described by an effective nonlinear 
Schrodinger (NLS) equation, was reported. Notice that in 
earlier studies on RH-TL models it was shown that two (or 
more) solitons propagating with the same group velocity, can 
be described by a system of two (or more) NLS equations 
ifTTl (see also ifTSi for theoretical as well as experimental 
results). Such coupled NLS equations have been studied 
extensively in nonlinear optics and mathematical physics; 
see, e.g., Refs. |[T9ti2n and references therein. They are 
well-known to give rise to a variety of vector solitons, 
including bright-bright, bright-dark, and dark-dark ones. 

In this work, we study analytically and numerically the in- 
teraction between backward- and forward-propagating soli- 
tons in a nonlinear CRLH-TL. Our model is a nonlinear ver- 
sion of a generic CRLH-TL model (see, e.g., Refs. O (H): 
the considered nonlinear element in the unit cell of the TL 
is the shunt capacitor, which simulates the presence of a het- 
erostracture barrier varactor (HBV) diode 1 6 1 (the capacitance 
of the HBV diode depends on the applied voltage). Starting 
from the discrete lump element model of the CRLH-TL, we 
derive a nonlinear lattice equation. First, we study the linear 
regime and show that, for certain frequency bands, RH- and 
LH-modes can propagate with the same group velocity. Next, 
we treat the nonlinear lattice equation in the framework of the 
quasi-discrete (or quasi- continuum) approximation (see, e.g., 
|[T5l [22l |23l and 119 1 for a review): we thus seek for enve- 
lope soliton solutions of the nonlinear lattice model, character- 
ized by a discrete carrier and a continuum envelope and em- 
ploy an asymptotic multi- scale expansion method, to derive a 
system of two coupled NLS equations. Each of these equa- 
tions describes the evolution of the envelope of a backward- 
(LH-) and a forward-propagating (RH-) mode. A systematic 
analysis of the system of the NLS equations reveals the ex- 
istence -in certain frequency bands- of three different types 
of vector solitons: (a) a backward-propagating bright soli- 
ton coupled with a forward-propagating bright soliton, (b) a 
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backward-propagating bright soliton coupled with a forward- 
propagating dark soHton. (c) a backward-propagating dark 
soHton coupled with a forward-propagating bright soliton, and 

The above analytical predictions are then tested against di- 
rect numerical simulations, which are performed in the frame- 
work of the original nonlinear lattice model. The results of the 
simulations verify the existence of the aforementioned types 
of vector solitons in the full TL model, but also offer important 
information regarding their robustness. In particular, results 
of direct simulations performed for long times indicate that 
bright-bright solitons are the most robust among the members 
of the vector soliton family. Indeed, the mixed (dark-bright 
or bright-dark) types are found to be less robust; however, 
the dark-bright solitons in a specific frequency band, although 
they are deformed during their evolution, are found to be more 
robust than those in other bands, as well as the bright-dark 
solitons, which are destroyed for the same propagation time. 
In any case, our results indicate the existence of all three types, 
robustness of bright-bright solitons and partial or substantial 
deformation of the other types. We can thus conclude that 
bright-bright (LH-RH), as well as dark-bright (LH-RH) soli- 
tons in certain frequency bands, have a better chance to be 
observed in experiments. 

The paper is organized as follows. In Section II, we intro- 
duce the nonlinear CRLH-TL model and the pertinent lattice 
equation, and derive the system of the two coupled NLS equa- 
tions (relevant details are also appended in an Appendix). In 
Section III, we present analytical and numerical results for 
each type of vector soliton. Finally, in Section IV, we summa- 
rize our conclusions. 



II. THE MODEL AND ITS ANALYTICAL 
CONSIDERATIONS 

A. The nonlinear CRLH-TL model 

We consider a generic CRLH-TL, composed by both right- 
and left-handed elements, as shown in its unit-cell circuit 
shown in Fig. [T] O O. The (RH) elements of this TL are 
the inductance Lr and capacitance Cr, while the LH ones 
are the inductance Ll and capacitance Cl- We assume that 
the TL is loaded with a nonlinear capacitance {Cr, while the 
capacitance Cl will be assumed to be fixed and voltage in- 
dependent). This can be implemented by proper insertion of 
diodes in the TL (see, e.g., pertinent experiments as well as 
theoretical work in Refs. CoHEl); in other words, we assume 
that the shunt capacitor Cr is nonlinear (see details below). 

Let us now consider Kirchhoff 's voltage and current laws 
for the unit-cell circuit of Fig.[T] which respectively read: 



Vn-l = Vn^LR^^Un 

at 



(1) 



In = In+l^lL^J^{CRVn), (2) 

where Un is the voltage across the capacitance Cl and II 
is the current across the inductor L^. The above equations. 
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FIG. 1 : The unit-cell circuit of the nonlinear CRLH model. 



together with the auxiliary equations Vn = LLdlL/dt and 
In = CLdUn/dt, lead to the following system: 

LrLlCl ^ {CRVn) ^Ll^ {CRVn) + LrCl ^ 



-LLCL^iVn+l + Vn-l - 2^) + K = 0. 



(3) 



To proceed further, we now consider a specific voltage- 
dependence for the nonlinear capacitance Cr. Here, we will 
assume that - for sufficiently small values of the voltage Vn 
- the function CR{Vn) can be approximated as follows, via a 
Taylor expansion: 

CRiVn) ^ Cro + CRoiVn - Vo) + IcRoiVn - Fo)^ (4) 

where Cro = Cr{Vo) is a constant capacitance correspond- 
ing to the bias voltage Vb, while C^q and C^q also assume 
constant values, depending on the particular form of Cr{V). 
Below, we will further discuss this approximation, in connec- 
tion with the HB V diode, used in the experiments described in 
Ref. 1 13 1 (similar varactor-type diodes were also used in the 
experiments of Ref. |7 |). 

Next, substituting Eq. ^ into Eq. ^ and using the scale 
transformations t ujsht [where uj^^ = {LlCro)~^] and 
K ^ [C^o(2C'i?o)"^]K, we obtain: 

r(K+i+K-i-2K) + (l + ^') 



d'^Vn 2 

dt^ ^ dt^' 

^2t/-2 j2t/-3 ^4t/-2 



d'Vn 



dt^ 



d^V^ 

^ ^ n 



df^ 



df^ 



(5) 



where the constant parameters 5, (3 and ^ are given by: 



5=^, p = lfL, ^ = ?^Cho. (6) 

Jsh Jsh oLy 



In the above expressions, /se and /sh denote series and shunt 
frequencies, while Jrh denotes the characteristic frequency 
related to the RH part of the unit-cell circuit, respectively; the 
above frequencies are defined as: 



/se — 



1 



2lly/LRCL 
1 



7 /sh — 



2t:^LlCrq ' 



7rh 



2lly/LRCR0 



(7) 
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Note that if /se//sh = 1, i.e., ^ = 1, then the CRLH-TL 
is usually referred to as balanced, in the sense that the char- 
acteristic impedances of the purely LH- and RH-TL, defined 
as Zl = and Zr = ^JLrICr^, are equal, i.e., 

Zl = Zr |2|. On the other hand, if f^e/ fsh > 1, i.e., ^ > 1, 
the LH part of the TL dominates, in the sense that the TL has 
a more pronounced LH behaviour (the serial branch features a 
capacitive character while the shunt branch an inductive one). 
In the opposite case, /se//sh < 1, i.e., S < 1, the RH part 
of the TL dominates and the TL has a more pronounced RH 
behaviour (the serial branch features an inductive character 
while the shunt branch a capacitive one). 

It is now relevant to adopt physically relevant parameter 
values for Eq. ([5|. For applications in the microwave fre- 
quency range (e.g., for microstrip lines ||2J or coplanar waveg- 
uide structures loaded with SRRs (31 - cf. also Ref. |23 1 for 
recent work), typical values of the capacitances and induc- 
tances involved in the CRLH structure are of the order of pF 
and nH, respectively. Here, we will use the values Lr = 1 nH, 
Cl = 0.1 pF, and Ll = 0.12 nH; thus, the frequencies in 
Eqs. take the values f^e = 15.92 GHz, /^h = 14.53 GHz 
and /rh = 5.03 GHz. On the other hand, as concerns the pa- 
rameters involved with the nonlinear capacitor Cr, we assume 
that the pertinent capacitance corresponds to a HBV diode, 
which is characterized by the following equation (Sj (see also 
ifTl . where the same form of C{V) is used, but different pa- 
rameter values): 



8=1.0954 



C{V) = CjoAda 1 



(8) 



where Cjo = 1.53 fF//im^ is the capacitance corresponding 
to bias voltage Vq = 0.2 V, A^a = 650 /im^ is the device 
area, V^r = 12 V is the breakdown potential, and the exponent 
m = 2.7 results from fitting experimental data. It is clear that, 
for sufficiently small V, by Taylor expanding Eq. ^ one ob- 
tains Eq. (|4]), where the constant parameter values involved are 
Cro = 1 pF, Crq = -0.24 pF/V and Crq = -0.08 pFA^^ 
To this end, the values of the normalized parameters S, (3 and 
fi appearing in Eq. ^ take the following values: 



1.1, /3^0.35, /i: 



-0.9. 



(9) 



Below, we will use these values for the purposes of our an- 
alytical and numerical considerations (we have checked that 
other values lead to qualitatively similar results). Notice that 
our choice leads to (5 > 1, i.e., we consider the case where the 
TL has a more pronounced LH character; however, when con- 
sidering the linear setting (see next subsection), this parameter 
will also assume other values, corresponding to the balanced 
and RH-dominated behaviour as well. 



B. Linear analysis 

We now assume plane wave solutions of Eq. ([5]), of the form 
Vn = Voex.-p[i{kn — ut)], where k and u denote the wave 
number and angular frequency, respectively, while the ampli- 
tude of the wave is 1. Substituting the above ansatz into 
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FIG. 2: (Color online) The dispersion relation showing the normal- 
ized frequency / / /sh as a function of the wave number k (in rad/cell) 
for different values of 3, i.e., 5 = 1.0954 (top panel), 5 = 1 
(middle panel), and 5 = 0.7746 (bottom panel). The solid (red) 
and dashed (blue) lines show the dispersion relation in the LH- and 
RH-frequency regions, respectively; RH± and LHi denote branches 
with A:>Oor/c<0. If^T^ la gap is formed; the width of the gap 
is |5 — 1| for 5 > 1 (top panel) ox 5 <1 (bottom panel). 



Eq. ([5]), and keeping only the linear terms in Vb, we obtain the 
following linear dispersion relation: 



4/3^ sin^ 



0. (10) 



The above result is illustrated in Fig.|2j where we plot the fre- 
quency ///sh as a function of the wave number k (in rad/cell), 
for three different values of 5 (note that here we consider one 
period of /c, i.e., —it < kj < tt). It is clear that for S — 1.0954 
(top panel) there exist two frequency bands where EM wave 
propagation is possible: the RH-band [high-frequency band 
depicted by dashed (blue) line], for 1.0954 < / < 1.4535, 
and the LH-band [low-frequency band depicted by solid (red) 
line], for 0.7538 < / < 1. In the same case = 1.0954), 
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there exists a gap for 1 < // f^h < where EM wave propa- 
gation is not possible. 

In the case where 6 = 1 (corresponding, e.g., to the value 
Cl = 0.12 pF) the gap vanishes (cf. middle panel of Fig. [2]) 
and the TL is balanced. In the balanced case, EM wave prop- 
agation is possible in two frequency bands as well: the RH- 
band [high-frequency band - cf. dashed (blue) line] with 
1 < ///sh < 1.405 and the LH-band [low-frequency band 
- cf. solid (red) line] with 0.7117 < ///sh < 1. 

Finally, for 6 = 0.7746 (corresponding, e.g., to Cl = 
0.2 pF), a gap appears again for ^ < ///sh < 1 (bottom panel 
of Fig. [2]). In this case too, there exist a RH-frequency band 
and a LH-frequency band, for 0.588 < ///sh < 0.7746 and 
1 < ///sh < 1.317, respectively. Note that in all cases, the 
RH± and LH± branches correspond to positive or negative k, 
respectively. 

Thus, generally, in the linear setting - and for a given fre- 
quency - the EM waves may either propagate in the RH region 
(forward wave propagation) or in the LH region (backward 
wave propagation). However, in the nonlinear setting, cou- 
pling between modes propagating in the LH and RH regime is 
possible (see, e.g., relevant earlier work in Refs. (TTIITJI). Be- 
low we will demonstrate that this is the case indeed, and study 
the coupling (interaction) between LH and RH modes with 
equal group velocities. Since the latter are tangents in the dis- 
persion curves, inspection of Fig. [2] shows that it is possible to 
identify domains, belonging to the RH± and LH^ branches, 
exhibiting parallel tangents, i.e., equal group velocities. 



8=1.0954 



To further elaborate on this, we may use Eq. ([TO]) to obtain 
the group velocity Vg = duj/dk: 



uj'^P^ sin/c 



(11) 



In Fig.|3] we show the dependence of the group velocity Vg on 
the normalized frequency f / fsh, for the values of 6 used in 
Fig. [2] Notice that the figure depicts only the group- velocity 
branches with Vg > 0- see solid (red) and dashed (blue) lines 
- corresponding, respectively, to the LH_ and RH+ branches 
of the dispersion curves; the branches with Vg < (pertinent 
to the LH+ and RH_ branches of the dispersion curve) are 
mirror symmetric with respect to the ones shown in the figure, 
due to the parity of the dispersion relation. 

Considering a horizontal cut of the group- velocity curves, 
say at Vg = 0.1 or Vg = 0.075 (see horizontal lines in the 
top and bottom panels of Fig. [3]), it is readily observed that, 
indeed, a LH_ and a RH+ mode can share a common group 
velocity (and interact in the nonlinear regime, as mentioned 
above). In fact, inspection of the group- velocity curves, say 
in the top panel of Fig. [3j shows that the maximum possible 
common Vg is given by Vg^^^ = 0.1339, the local maximum 
of Vg, occurring at / = 0.9391, in the (shorter in height) LH 
frequency band. Then, one can divide each of the LH and 
RH group- velocity curves into two sub-regions, depending 
on the sign of the group- velocity dispersion (GVD), dvg/duj, 
where such coupling with equal group velocities may occur. 
These subregions are: (a) the sub-bands I (0.7538 < ///sh < 
0.9391) and II (0.9391 < ///sh < 1) for the LH-frequency 
band, characterized by positive and negative GVD respec- 
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FIG. 3: (Color online) The group velocity Vg as a function of the 
normalized frequency ///sh (for 6 = 1.0954). The solid (red) and 
dashed (blue) lines indicate branches corresponding to the LH_ and 
RH+ regimes, respectively. The intersection of the group velocity 
curves with the horizontal solid (black) line depicts frequencies of 
modes with the same group velocity, Vg = 0.1. Regions I, II, III and 
IV indicate possible interactions between LH_ and RH+ modes with 
the same Vg but different signs of GVD. 



tively, and (b) the sub-bands III (1.0954 < ///^h < 1.1195) 
and IV (1.356 < ///sh < 1.4535) for the RH-frequency 
band, again characterized by positive and negative GVD re- 
spectively. Thus, nonlinear LH and RH modes of equal Vg 
can feature the following four different possible interactions: 

1 . LH-mode in band II and RH-mode in band IV, both fea- 
turing negative GVD. 

2. LH-mode in band I and RH-mode in band IV; here, the 
LH (RH) mode features positive (negative) GVD. 

3. LH-mode in band I and RH-mode in band III, both fea- 
turing positive GVD. 

4. LH-mode in band II and RH-mode in band III; here, the 



5 



LH (RH) mode features negative (positive) GVD. 

It is clear that the above set of possibilities arises from the 
existence of the gap in the considered case with S = 1.0954. 
A similar situation also occurs for S < 1, e.g., for S = 0.7746 
as in the bottom panels of Figs.|2]and[3] On the other hand, for 
(5 = 1 the gap does not longer exist and, thus, the only possi- 
ble interaction is between a LH-mode with positive GVD and 
a RH-mode with negative GVD; this interaction can occur for 
group velocities Vg < 0.175, i.e., beneath the dashed horizon- 
tal line in the middle panel of Fig. [3] This possibility, however, 
is already taken into regard - cf. case (2) above; furthermore, 
soliton formation in the balanced CRLH-TL (5 = 1) has al- 
ready been studied in the literature |11|. For these reasons, 
below we will proceed by analyzing the case corresponding 
to S = 1.0954, which offers all possible scenarios; it is clear 
that the case of ^ = 0.7746 shares similar qualitative features; 
this similarity extends beyond the linear wave case and into 
the nonlinear solitonic one. 

Although, as explained above, we are not going to analyze 
soliton formation and soliton in the special case of the bal- 
anced CRLH-TL with 6 = 1, it is worth mentioning the fol- 
lowing. In the case of (5 = 1, the dispersion relation exhibits a 
Dirac point, namely it is approximately linear in the vicinity of 
k = 0, i.e., ^ ±[1 + {f3/2)k] - cf. middle panel of Fig.[2] 
The emergence of Dirac points is particularly interesting in 
the two-dimensional (2D) setting of triangular and hexago- 
nal lattices arising in different contexts, such as optics f24\ , 
atomic Bose-Einstein condensates | 25 1, and the so-called pho- 
tonic graphene 1261 . This has also led to an interest in this 
subject from a rigorous mathematical perspective |27|. It is 
thus quite interesting that, in principle, 2D balanced CRLH- 
TLs may host a variety of fundamental effects, such as conical 
diffraction, formation of topological defects, and even phase 
transitions, as in Refs. I.24.-.26 J . 



(where Vg is the common group velocity) and T = e'^t, while 
exp(z6>j), with Oj = kjU — Ujt, are the (discrete) carriers 
of frequencies ujj and wavenumbers kj . Finally, e is a formal 
small parameter setting the field amplitude and the slow scales 
of the envelope functions. 

At this point, we should note that the field Vn as expressed 
in Eq. ([12]) is, in fact, the leading-order form of a more gen- 
eral ansatz employing multiple time and space scales. In this 
context, use of a formal multi- scale expansion method leads 
to a hierarchy of equations at various powers of e, which are 
solved up to the third-order. Here, we will present the main 
results and provide further details in the Appendix A. Particu- 
larly, from the first- and second-order problems [i.e., at orders 
0{e) (linear limit) and O(e^), respectively] we derive the dis- 
persion relation, Eq. ( [TO] ), and the group velocity, Eq. ( [TT] ). 
Finally, at the next order, (9(e^), we obtain the following cou- 
pled NLS equations: 

idrV, + ^D^dj^V, + (^iilFip +^i2|F2p) Vi = 0, (13) 

idTV2 + lD2dj,V2 + (^2l|^l|' +^22|^2|') ^2 = 0, (14) 

where the normaUzed GVD coefficients Dj, the self-phase 
modulation (SPM) coefficients gjj, and the cross-phase mod- 
ulation (CPM) coefficients gj,3-j (with j = 1, 2) are respec- 
tively given by: 



933 



93,3-3 



cot ki 



5^) 



(3/x - A^) , 



2(^4 _ ^2) 



(6m - Ba-i) , 



(15) 
(16) 
(17) 



C. Nonlinear analysis: the coupled NLS equations 

To describe the coupling between a RH and a LH nonlin- 
ear mode with equal group velocities, we will use the quasi- 
discrete approximation, which takes into regard the inherent 
discreteness of the system (see, e.g., Ref. 1 19] for a review, 
and Refs. |15, 23] for relevant recent work). Generally, this 
approach allows for the description of quasi-discrete envelope 
solitons (usually satisfying an effective NLS model), charac- 
terized by a discrete carrier and a slowly-varying continuum 
envelope. In our case, since we are interested in the descrip- 
tion of two different modes, we seek for a solution of Eq. ([5]) 
in the form: 



2 

^ ^jn{X, T) exp{iOj) + c.c, 



(12) 



In Eq. (12) 



denotes complex conjugate 
= 1,2 correspond to the LH and RH mode 



where "c.c' 
subscripts j 

Vjn{X^r) are unknown (continuous) slowly- varying enve- 
lope functions depending on the slow scales X = e{n — Vgt) 



and the coefficients Aj and Bs-j are defined in Appendix 
A. Next, using scale transformations, we measure normalized 
time T and densities |V^p in units of and \Di/gjj\ 

respectively, and cast Eqs. ([T3])-([T4]) in the form: 



idrVi 
idTV2 



djcV2 ^ {X2\Vi\' 



Xi\V2nVi 



■(T2\V2nV2 



0, (18) 
0, (19) 



where 



s = sign(i:>i), cr^- = sign (^^., ), 

D2 , 912 . 921 



d 



Ai 



1^221 



1^11 1 



(20) 



As seen from Eqs. ( [18] ), in the absence of coupling (Aj = 0) 
the evolution of either the LH mode Vi or the RH mode V2 
is described by a single NLS equation; the latter, supports 
soliton solutions of the dark or the bright type, depending on 
the relative signs of dispersion and nonlinearity coefficients 
(see, e.g., Ref. |20|). Particularly, the mode Vi (F2) supports 
dark solitons for scri < (da2 < 0) or bright solitons for 
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sdi > (da2 > 0). These conditions, however, are modified 
for Xj 7^ and various types of coupled (aUas vector) soH- 



tons can be found in the full version of Eqs. ( 18 ). Below we 
will present these types of coupled backward- and forward- 
propagating solitons, belonging, respectively, to the LH and 
RH frequency bands. 

Before proceeding with the presentation of the coupled soli- 
ton solutions we make the following comments. First, in some 
cases, solitons will be found in a stationary form. However, 
using these stationary solutions, one can also find travelling 
soliton solutions, with an additional free parameter, i.e., the 
velocity C, by means of the following Galilean boost: 



Vi{X,T) Vi{X-CT,T) 



X exp 



{i 



cx 



V2{X,T) V2{X-CT,T) 



X exp • 



d 



CX 



C2 



']} 



(21) 



(22) 



Second, it is interesting to note that, contrary to what is often 
the case in the mathematically studied multi-component vari- 
ants of the NLS equation |21 1, the model of Eqs. ([T8])-([T9) 
does not necessarily respect the condition Ai = A2. The latter 
condition ensures the existence of an underlying Hamiltonian 
structure and is customary in other physical applications (such 
as atomic physics |28|). Nevertheless, as we will see below, 
this is not a necessary condition for the existence of the exact 
solutions considered below. 



equal to 0.01, we have fixed the value of the small parameter 
as e = 0.02, and we have used periodic boundary conditions. 
Use of the latter leads to the requirement that the wavenumber 
/c of a dark soliton component must be equal to 27rq/p, with 
q^p e Z and q also being odd. 

In all figures below (Figs.[5p5]), except if stated otherwise, 
we show the density plots of Vn, the spatial profile of Vn at 
t = 2000, as well as the time evolution of the center of mass 
X{t) and the quantity P{t). 

Regarding the evolution time of the simulations, we should 
note the following. Most of our simulations are performed for 
relatively large normalized times - typically up to t ~ lO*" 
in some cases. However, given our time normalization, the 
physical unit time (set by the frequency fsh = 14.529 GHz) 
is very small, namely to = {27r fsh)~^ ^ H picoseconds 
(see Sec. II. A). Actually, since all characteristic frequencies 
of the system (see Eq. ([7])) are in the microwave regime, all 
characteristic times are less than a nanosecond and, thus, ob- 
viously, simulations for time t even of the order of 10^ are 
extremely time-consuming. Nevertheless, our results for nor- 
malized times up to t = 10^ (corresponding to a physical time 
of the order of a tenth of millisecond), demonstrate a good 
agreement with our analytical predictions in suitable cases 
(see below). Furthermore, the results of such long simula- 
tions can also be used as a reliable indication of the solitons' 
robustness. Hence, in the case where the solitary waves are 
found to be very robust, we expect that they would survive for 
the longer time scales that would render them experimentally 
observable. 



III. SOLITON INTERACTIONS IN DIFFERENT 
FREQUENCY BANDS. NUMERICAL RESULTS 

A. Numerical procedure. 

Let us now proceed to study numerically the evolution of 
the coupled solitons presented in the previous section in the 
framework of the fully discrete model of Eq. ([5]). 

In order to compare the analytical approximations with the 
results of numerical simulations, we will make use of two di- 
agnostic quantities: the first one is the evolution of the center 
of mass defined as: 



X{t) = 



Z^f} — - '^^r, 



En 
n 



-N ' 
n=N 
-N 



1/2 



(23) 



and the second one, is a power-like quantity defined as: 

n=N 

m= E (24) 

with 2N + 1 being the lattice size. The above quantities can 
readily be determined for each type of vector solitons that is 
predicted analytically in the framework of the coupled NLS 
equations. 

In all simulations, which have been performed by means of 
a fixed-step 4th-order Runge-Kutta scheme with a time step 



B. Bright-bright solitons in bands II and IV. 

First, we consider the interaction between a backward prop- 
agating soliton, with a frequency lying in band II, and a for- 
ward propagating soliton, with a frequency lying in band IV. 
In this case, s = —1 (cf. Fig. [Sj, while the other dispersion 
and nonlinearity coefficients are shown in Fig.|4]as functions 
of the normalized frequency f / fsh (for S = 1.0954). It is ob- 
served that (Ji = (72 = — 1, and also Xi ^ X2 = A, while 
the dispersion coefficient d takes values d < — 0.25. Thus, to 



a first approximation, the system of Eqs. ( 18 )-( 19 ) takes the 
form: 

idTVi-^dj,Vi^{X\V2\^ -\Vi\^)Vi=0, (25) 

idTV2^^dj,V2^{X\Vi\^ -\V2\^)V2 = 0. (26) 

where A < and d < 0, as can be seen in the top panel 
of Fig. |4] The above system is generally non-integrable 
and soliton solutions can not be found in an explicit analyt- 
ical form. However, there exists a specific frequency value, 
namely f / fsh — 0.9654 where the dispersion coefficient d 
take the value d ^ —1 and the nonlinearity coefficients A 1,2 
take the values Ai = A2 = —1.7 (see the intersection point 
of the relevant curves depicted by a star in the top panel of 
Fig. |4]). This case corresponds to a (common for both com- 
ponents) group velocity Vg = 0.1288, which occurs when 
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FIG. 4: (Color online) Parameters for soliton interactions in bands II 
and IV. Top panel: the dependence of the parameters Ai [thin solid 
(blue) line], A2 [dashed (blue) line] and d [bold solid (red) line] on 
the normalized frequency ///sh- Bottom panel: the nonlinearity co- 
efficients gii [solid (red) line] and ^22 [dashed (blue) line] as func- 
tions of ///sh- The parameter 6 takes the value S = 1.0954. Stars (in 
black) in both panels corresponding to ///sh = 0.9654, for which 
d — (71,2 = — 1 and Ai = A2 = A = — 1.7. These values are used 
for the simulations shown in Figs.|5]and[6]below. 



the (normalized) carrier frequencies for the modes Vi and V2 
take, respectively, the values /i//sh = 0.9654 (as mentioned 
above) and /2//sh = 1.3653. Then, symmetric bright-bright 
standing soliton solutions can be found in the following form 
(see, e.g., Ref. |29|): 



Vi=V2 



2i 



sech{V2£X) exp(-i£T), (27) 



where £ is an arbitrary parameter. Using the above expres- 
sions, we can now approximate the unknown voltage Vn{t) in 
Eq. ([5]), in terms of the original coordinates n and t, as fol- 
lows: 



Vn{t) ^ Vo[Ri{n,t)cos{kin-nit) 
+ R2{n,t) cos{k2n — Q2t), 

where functions Ri and R2 have the following form: 



Ri =sech[V2£e{n-Vgt)], 
R2 = \ l ^ sech[V^e(n - Vgt)], 



922 



(28) 

(29) 
(30) 



while the solution amplitude Vq and the frequencies Qj (j = 
1, 2) are given by: 



Vn 




(31) 
(32) 



with Uj = fj/fsh- Now, substituting Eq. (28) into Eqs. (23) 



and ( 24 ) and supposing that e is small enough, we obtain for 
our diagnostic quantities: 



X{t) 
P(t) 



Vgt, 

1/2 



922 



(33) 
(34) 



In Figs. [5] and |6] we show the outcome of the simulations 
for short and long times, respectively, of a bright-bright soli- 
ton with £ = 1 and N = 500. The parameters used are 
/i = 0.96545 and /s = 1.36535, which gives h = -0.4061 
and k2 = 1.8576, i.e., a bright-bright soliton in bands II and 
IV. In Fig. [5] it is evident that the agreement between analyti- 
cal and numerical results pertaining to the soliton profile, and 
the evolution of the center of mass and power diagnostics, is 
very good. In the case shown in Fig. [6j we have performed 
a very long simulation, up to normalized times t = 10^. It 
is clear that that the initial pulse does not spread out, which 
indicates the soliton robustness: the top panels of the figure - 
and particularly the snapshots of the pulse profile at t = 10^ 
- clearly show that the soliton persists as a stable object up to 
the end of this long simulation time. We note in passing that a 
fragment of the soliton is backscattered when the soliton starts 
its motion at t = (due to the approximate nature of our ana- 
lytical solution profile). Notice that despite this emission and 
the subsequent interaction of the fragment with the "distilled" 
solitary wave, the coherent structure remains robust and pre- 
serves its characteristics throughout the evolution thereafter. 



C. Bright-dark solitons in bands I and IV 

Next, we consider the interaction between a backward prop- 
agating soliton, with a frequency lying in band I, and a for- 
ward propagating soliton, with a frequency lying in band 
IV. In Fig. [7] we show the dependence of the parameters Ai, 
A 2 and d (top panel), as well as of the nonlinearity coef- 
ficients (bottom panel), on the normalized frequency ///sh 
(for S = 1.0954). In this case, 5 = +1 (cf. Fig.[^, while 
(Ji = a2 = —I and, thus, Eqs. ([T8])-([T9]) are reduced to the 
form: 



idrVi + -dj,Vi + (A1IF2I' - l^iP) = 0, (35) 

idTV2 + ^dj,V2 + - IF2I') V2 = 0, (36) 

where A 1,2 < and d < 0, as can be seen in the top panel 
of Fig.|7] The above equations are no longer of the Manakov 
type and, thus, generally, they are not completely integrable. 
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FIG. 5: (Color online) Bright-bright solitons in regions II-IV. Top 
left: density plot of the space-time evolution of Vn obtained numer- 
ically. The top right panel compares the analytical and numerical 
profiles of at t = 2000. The bottom panels show the time evolu- 
tion of the center of mass (left) and the power diagnostic (right). The 
parameters used are /i = 0.96545 and /2 = 1.36535, which gives 
ki — —0.4061 and /c2 = 1.8576, i.e. a bright-bright soliton in bands 
II and IV (this particular choice corresponds to the points depicted by 
stars in Fig.|4]). The difference in the powers can be attributed to the 
approximate nature of our solution. 



S=-l-1 

Q — 




FIG. 7: (Color online) Same as Fig.|4] but for soliton interactions in 
bands I and IV. Stars depict parameter values used for the simulations 
shown in Figs.[8]and[9]below. 



propagating dark soliton, whose exact analytical form is: 





FIG. 6: (Color online) The bright-bright soliton of Fig.^is evolved 
until t — 10^. All the panels are similar to that of Fig. ^except for 
the top right. In this panel, snapshots of the soliton at t = 2 x 10^ 
and t = 10^ are compared to the initial condition of the simulation 
in order to examine its robustness under a very long evolution time. 



Nevertheless, standing wave solutions can still be found in the 
form of coupled bright-dark solitons, with the frequency of 
the bright (dark) soliton component being in the LH (RH) fre- 
quency band. Therefore, here we have a case of coupled soli- 
tons, a backward-propagating bright soliton and a forward- 



Vi(X,T) 



\d^\2\) 
/7^tanh(6X) exp(— zi^2^), 



sech(6X) exp(— zi^iT), (37) 
(38) 



where the soliton amplitude parameters vi^2 and the inverse 
width h are connected via the following equations: 



z/i = - 



^2 



2(d + A2) 
^2 



[l + (2c^ + A2)Ai], (39) 
(I-A1A2), (40) 



with 1 — A1A2 < in the considered LH frequency band (no- 
tice that (i + A2 < as well). It is thus clear that the above 
solutions are characterized by one free parameter. 

Employing the solutions ([37]) -([38]), we can again approxi- 
mate a solution of Eq. (|5| for the voltage Vn(t), in terms of 
the original coordinates n and as follows: 



Vn{t) ^ Vo[^i{n,t)cos{kin-nit) 
+ '^2{ri,t)cos{k2n - Q2t)], 



(41) 



where 



^1 = ^1,0 sech[e6(n - Vgt)], ^1,0 = 
^2 = ^2,otanh[e6(n - Vgt)], ^2,0 
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FIG. 8: (Color online) Bright-dark solitons in regions I-IV. Top left: 
density plot of the time evolution of Vn obtained numerically. The 
top right panel compares the analytical and numerical profiles of Vn 
at t = 2000. The bottom panels show the time evolution of the center 
of mass (left) and the power diagnostics (right). The parameters used 
are fi = 0.8831 and = 57r/8 ^ 1.9625, which gives ki = 
— 1.0404 and /2 = 1.3748 (see corresponding points depicted by 
stars in Fig. [7]), i.e. a bright-dark soliton in the I and IV bands. 
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FIG. 9: (Color online) The bright-dark soliton of Fig.j^is evolved 
until t = 2 X 10^. All the panels are similar to that of Fig^except for 
the top right. In this panel, snapshots of the soliton at t = 2x10^ and 
t = 5 X 10^ are compared to the initial condition of the simulation. 
Notice that the center of mass is not bounded into [—N, N] ; this can 
be caused by the soliton splitting. 



In this case, the solution amplitude Vb and the frequencies Qj 
(j = 1,2) are given by: 



Vn 



911 



Qj = ujj^e^Vj\Di\. 



(44) 
(45) 



In order to get an expression for the center of mass similar 



to that of the bright-bright soliton ( 34 ), we must define it as: 



Xd{t) 
where 



- ^*^,o cot(fc2) sin(2f22t), (46) 



-2,0 



eh 



m = Ko + Koi^^^ - 1) + ^^2,0[l - COs(2^]2t)]. 

(47) 

Substituting Eq. ( [4T] ) into Eqs. ( [46| ) and ([24]) we can once 
again obtain relevant expressions (provided that e is small 
enough) for the center of mass and power: 



1/2 



(48) 



(49) 



Figures [8] and [9] show the evolution of a bright-dark soliton 
(and its characteristics) in bands I and IV with U2 — ^ and 
N = 1220. The parameters used are /i = 0.8831 and k2 = 
57r/8 ^ 1.9625, which give ki = -1.0404 and /s = 1.3748. 
In this case, it is clear that although bright-dark solitons do 
exist, the agreement between analytical and numerical results 
becomes worse over time. Also, as shown in the top right 
panel of Fig.|9] the pulse profile indicates that the bright-dark 
soliton is not a robust object. In particular, at time t = 5 x 10^, 
it is clear that the configuration has dramatically changed its 
character. 



D. Dark-bright solitons in bands I-III and II-III 

Finally, we consider the cases of coupled solitons in bands 
I and III, and also in bands II and III. In both cases, as is 



observed in the top panels of Figs. 10 (bands I-III) and 11 
(bands II-III), we have that the parameter A2 <C 1. Now, the 
NLS equations for Vi and V2 take the following form: 

idrVi + ^dj,Vi + (A1IF2I' - l^il') Vi = 0, (50) 

idrVi + ^dj,Vi + - V2 = 0, (51) 

where s = +l(s = —1) corresponds to solitons in bands I 
and III (II and III). Note that, in this case, ai = (72 = — 1, 



while d > and Ai,2 < 0, as shown in Figs. 10 and 11 As in 



the previous case, the equations (50)-(51 ), generally, are not 
completely integrable. Nevertheless, standing wave solutions 
can still be found in the form of coupled dark-bright solitons, 
with the frequency of the dark (bright) soliton component be- 
ing in the LH (RH) frequency band; therefore, here we have 
a case of a backward-propagating dark soliton, coupled with 
a forward-propagating bright soliton, whose exact analytical 
forms are: 



Vi(X,T) = v^tanh(BX)exp(-%T), 



(52) 



V2{X,T) = j!?lil±^sech(BX)exp(-z7y2T),(53) 
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FIG. 10: (Color online) Same as Fig.|4] but for soliton interactions in 
bands I and III. Stars depict parameter values used for the simulations 
shown in Figs. 1 1 2|and| 1 3 |below. 
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FIG. 11: (Color online) Same as Fig.|4] but for soliton interactions 
in bands II and III. Stars depict parameter values used for the simu- 
lations shown in Fig s . \i4\ and p3] below. 
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FIG. 12: (Color online) Dark-bright solitons in regions I-III. Top 
left: density plot of the time evolution of Vn obtained numerically. 
The top right panel compares the analytical and numerical profiles 
of at t = 2000. The bottom panels show the time evolution of 
the center of mass (left) and the width diagnostic (right). Parameters 
used are ki = — Gtt/S ^ —1.884 and /2 = 1.1002, which gives 
/i = 0.8003 and k2 = 0.1232 (cf. points depicted by stars in 
Fig.ITol, i.e. a dark-bright soliton in bands I and III. 
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FIG. 13: (Color online) The dark-bright soliton of Fig.[T2]is evolved 
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until t — 2 X 10^. All the panels are similar to that of Fig 
except for the top right. In this panel, snapshots of the soliton at 
t = 2xl0^andt = 5xl0^ are compared to the initial condition of 
the simulation in order to examine its robustness under a very lengthy 
time evolution. 



where the soliton amplitude parameters 771^2 and the inverse 
width B are connected via the following equations: 



2{dXi^s) 



[d+(dAi+25)A2], (54) 
(I-A1A2). (55) 



Following our previous considerations, we may again use the 
above solutions and approximate the voltage Fn(^) in Eq. ([5|, 
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FIG. 14: (Color online) Dark-bright solitons in regions II-IIL Top 
left: density plot of the time evolution of Vn obtained numerically. 
The top right panel compares the analytical and numerical profiles 
of at t = 2000. The bottom panels show the time evolution of 
the center of mass (left) and the width diagnostic (right). Parameters 
used are ki = — 37r/23 ^ —0.4095 and /2 = 1.1162, which gives 
/i = 0.965 and k2 = 0.2758 (cf. stars in Fig.[TTJ, i.e. a dark-bright 
soliton in bands II and III. 



in terms of the original coordinates, as follows: 

Vn{t) = Vo[^i{n,t)cos{kin-nit) 
+ ^2{n,t)cos{k2n - Q2t)], 



(56) 



where functions $i and ^2 are given by: 

= tSinh[eB{n - Vgt)], 
^2 = ^2,0 sech[eB{n - Vgt)], 



^2,0 



l{d^sX2) 
\dXi^s\ 



911 



922 



(58) 



while the rest of the soliton parameters are: 



Vo = 2eJr]i 



911 



(59) 
(60) 



Figures [12] and [13] show the outcome of the simulations for 
a dark-bright soliton in bands I and III (with r^i = 3 and N = 
3333), while Figs. [14] and [15] correspond to a dark-bright soli- 
ton in bands II and III (with r^i = 3 and N = 1553). The pa- 
rameters used are ki = — 67r/5 ^ —1.884 and /2 = 1.1002, 
which gives /i = 0.8003 and /c2 = 0.1232, in bands I and 
III, and ki = -Sir/ 23 ^ -0.4095 and = 1.1162, which 
gives /i = 0.965 and k2 = 0.2758, in bands II and III, re- 
spectively. In the latter case, the relatively large values of the 
number of particles and of r]i used are motivated by the neces- 
sity of a vanishing tail for the bright component at the edges 
of the lattice. As seen in this set of figures, dark-bright soli- 
tons in bands II-III and I-III do exist, as predicted in theory. 
Furthermore, it is observed that the former are less robust than 



FIG. 15: (Color online) The dark-bright soliton of Fig.[T4]is evolved 
2 X 10^ 



until t 



All the panels are similar to that of Fig. 



14 



[ except 

for the top right. In this panel, snapshots of the soliton at t = 10^ and 
t — 2 X 10^ are compared to the initial condition of the simulation. 
Notice that, as in Fig. [9] the center of mass is not bounded into 
[—N^ N] \ as in that case, the soliton splits at long evolution times. 



the latter, as seen both from their stronger deformation and the 
fact that they "lose" their solitary wave character earlier. This 
can be observed, e.g., in the strong fluctuations in the evolu- 
tion of the soliton center in the bottom left panel of Fig. 15 
or perhaps most notably in the substantial modification of the 
wave profile upon long propagation in the top right panel of 
the same figure. On the other hand, the dark-bright solitons of 
bands I-III seem to essentially preserve their structure even in 



(57) the long evolution of Fig. 13 



IV. CONCLUSIONS 

In conclusion, we have used both analytical and numeri- 
cal techniques to study the existence, stability and dynam- 
ics of coupled backward- and forward-propagating solitons 
in a composite right/left-handed (CRLH) nonlinear transmis- 
sion line (TL). The considered form of the TL was a quite 
generic one, finding applications to the modelling of a wide 
range of LH systems and devices, with "parasitic" RH behav- 
ior, such as resonators, antennas, directional couplers, among 
others |T-^|. 

Our analysis started with the derivation of a nonlinear lat- 
tice equation governing the voltage across the fundamental 
(unit cell) element of the transmission line. In the linear 
regime, we derived the dispersion relation for small- amplitude 
linear plane waves and showed that they may either propagate 
in a right-handed (RH) high-frequency region, or in a left- 
handed (LH) low-frequency region. We also identified fre- 
quency bands where RH- and LH-modes can propagate with 
the same group velocity. 

Using the above result, we then investigated the possibil- 
ity of nonlinearity-assisted coupling between LH- and RH- 
modes. This way, in order to analytically treat the nonlinear 
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lattice equation, we used the so-called quasi-discrete approx- 
imation. The latter is a variant of the multi-scale perturba- 
tion method, which takes into regard the discreteness of the 
system by considering the carrier (envelope) of the wave as 
a discrete (continuum) function of space. Employing this ap- 
proach, we derived, in the small- amplitude approximation and 
for certain space- and time-scales, a system of two coupled 
nonlinear Schrodinger (NLS) equations for the unknown volt- 
age envelope functions. This system was then used to predict 
the existence of coupled backward- and forward-propagating 
solitons, of the bright-bright, bright-dark and dark-bright type, 
respectively. 

The above existence results, as well as the propagation 
properties and the potential robustness of these vector soli- 
tons, were then investigated for each of the possible scenar- 
ios. This was done by means of direct numerical simulations 
of the full CRLH-TL nonlinear lattice model, using as initial 
conditions the analytical forms of solitons predicted by the 
perturbation theory. In the simulations, apart from the evo- 
lution of the shape, we also studied the evolution of the cen- 
ter of mass and a power-like quantity of the various solitons. 
Our numerical results have confirmed the existence of the var- 
ious types of solitons predicted analytically, but have also re- 
vealed their distinct robustness characteristics. In particular, 
we found that bright-bright solitons feature a robust propaga- 
tion over long times. On the other hand, as concerns solitons 
of the mixed-type (namely dark-bright and bright-dark ones), 
we found that, in specific frequency bands (bands I-III), dark- 
bright solitons are more robust than those in other bands (i.e., 
II-III) or bright-dark solitons: dark-bright solitons in bands 
II-III and bright dark solitons preserve their shape only for fi- 
nite times and, for sufficiently long evolutions, they are either 
destroyed (bright-dark) or are significantly deformed (dark- 
bright). 

We can thus postulate that from all types of solitons pre- 
dicted analytically, bright-bright and dark-bright ones (in 
bands I-III) are the most likely ones to be experimentally ob- 
servable. In all cases, our numerical results were found to 
corroborate the analytical predictions, at least up to the times 
during which the solitary waves propagate robustly. 

It would be interesting to study other types of nonlinear 
CRLH-TL lattice models modelling realistic structures com- 
posed by LH-metamaterials. In that regard, a pertinent in- 
teresting direction would be the investigation of the effects 
of damping and driving, which may lead to robust nonlin- 
ear waveforms which would constitute dynamical attractors in 
such settings. Additionally, the study of higher-dimensional 
settings is a particularly challenging problem. In the latter 
context, in addition to simpler (yet genuinely higher dimen- 
sional, or even quasi-one-dimensional) solitary wave struc- 
tures, more complex waveforms may be realizable such as 
vortices. The exploration of such states and their dynamical 
robustness will be reported in future publications. 
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Appendix A: The perturbation scheme 

Our analytical approximation relies on the use of the quasi- 
continuum approximation, which is a variant of the method of 
multiple scales |30|. We introduce new independent temporal 
variables, = (n = 0, 1, 2, • • • ), and accordingly expand 

the time derivative operator dt sls dt = dtg -\- edt^ -\- Next, 

we seek solutions of Eq. ^ in the form: 



c.c. 



(Al) 



£=1 



Then, we substitute Eq. ( Al ) into Eq. ([5| and employ a contin- 
uum approximation for the envelope functions Un, i.e., Un 
u{x), where x = na and a being the lattice spacing (the lat- 
ter parameter does not appear in the results below, as one may 
readily rescale x as x/a). Furthermore, we introduce the new 
spatial variables Xn = ^^x and, thus, dx = + edx^ + — 
To this end, equating coefficients of like powers of e, we ob- 
tain the following (first three) perturbation equations: 

0(e): Loiii=0, (A2) 
0{e^) : LoU2 + Lmi + NquI = 0, (A3) 
0{e^) : LiU2 + + No[uiU2 + /J^uf] = 0, (A4) 



where the operators are given by 



1 + + 4/3^ sin^ 



2 7 dtl 



(A5) 



dtldh 
2iP^ sin k 



21 



4/3^ sin^ 



2 J dtodti 



dt^dxi ' 



(A6) 



1 + (5^ + 4/3^ sin^ 



d^^^dtodt2 
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No 



dtldtl 
4z/3^ sin k 



2z/3^ sin k 



dtodtidxi 



dtldl 



dtl 



(Al) 
(AS) 



We now seek for a solution of the linear problem, Eq. ( A2 ), 
in the form: 



ui = ^Vj{xi,X2, . . . ,ti,t2, . . .) ex.p{iOj) +C.C., (A9) 

where subscripts j = I and j = 2 correspond to the LH 
and RH frequency bands, Vj is an unknown complex function, 
Oj = kjXQ — ujjto, while the wavenumbers kj and frequencies 



ujj satisfy the dispersion relation provided in Eq. (10). 
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Next, substituting Eq. ( |A9| ) into Eq. ( |A3| ), we obtain the 
non-secularity condition for / = 1: 



dVj 



sin kj 



2cj2- (1 + ^52+4/32 sin'^) 



^=0, (AlO) 
0x1 



which suggests that Vj = Vj{X,X2r ' ' 7^2,---)' where 
X = xi — Vg.ti, while the group velocitiesv^^ resuk self- 
consistently as Vg. — dujj / dkj [cf . Eq. (11)]. Employing 
Eq. (AlO), we may determine from Eq. (Aj), for / = 2, the 
unknown field U2 : 



U2 



E 



Gi{wi - uj2,ki - ^2) 



V,2exp(i26'^) 



-^Fj{xi,X2, - ■ ■ ,ti,t2,- ■ ■) + C.C. 

where functions Gj{ujj,kj) are given by: 



VtV2eMi{^i+^2) 
V^V^ eMi{Oi - 62) 

(All) 



G3 



(1 + (5^ + 4/32 sin2A;j-)(2wj-)^ 



X2 

{1 + 5^ + 4/3^ sm\'^ 



(A12) 



-))(W1+W2)2 



+ (Wi +W2)'' + (5^ 



(A13) 



G4 = 



(l + ,5"+4/?" sin^(- 
+ (wi - ^2)" + '5^- 



))(wi - W2)^ 



(A14) 



On the other hand, functions Fj{xi,X2, - ■ ■ , * i , ^2 , ■ ■ • ) can be 



derived at the order 0(e^), by means of the equation: 

i2M2 + N2UI = 0, 

which leads to the result: 



' 3 ~ , A 



(A15) 



(A16) 



To this end, we arrive at the following expression for U2 : 

2 



U2 



= -Yl ^^-^Z exp(i2^,) - CSV1V2 exp(z(^i + O2) 



where 



-04^1^2* exp(i(^i - ^2) - ^ co, |lS |2 + C.C, (A17) 

i=i 



C3 

C4 



4^2(4^2 _ ^2) 

Gj{2io„2k,) ' 

2[{0Ji+UJ2f-6^{0Jl+0J2f] 

G3(wi +W2,fci +^2) 
2[K-C^2)"-.52(c^l-C^2)'] 

G4(a;i - W2, fci - k2) 



_ 2c^2j2 



Finally, defining the coefficients: 



— Cqj + Cj , 
Bs-j = Co3-j + C3 + C4, 



(A18) 
(A19) 
(A20) 
(A21) 



(A22) 
(A23) 



and using the variables X = xi — Vgti = e{n — Vgt) and 
T = t2 = e^t, we derive from the non-secularity condition at 
0{e^) the coupled NLS equations (p^. 
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